R es un entorno y un lenguaje de programación estadística.
Es una de las herramientas más usadas en el análisis de datos. Es de código libre, su uso es gratuito y dispone de una comunidad muy activa que contribuye al desarrollo de herramientas y funciones adicionales (más de 10.000 paquetes) y a la resolución de preguntas a través de la red.
https://www.r-project.org/about.html https://www.rdocumentation.org/ https://stackoverflow.com/questions/tagged/r
R
Se puede instalar desde https://cran.r-project.org/.
RStudio
Es una interfaz gráfica (o IDE) para la realización de cálculos y la programación en R.
Se instala desde https://www.rstudio.com/, versión RStudio Desktop (Open Source License).
R puede usarse desde la consola introduciendo directamente las instrucciones correspondientes.
2 + 3 ## [1] 5
8 ^ 10## [1] 1073741824
log(100) # Logaritmo neperiano## [1] 4.60517
# indica que lo que sigue es un comentario
Un script en R es un conjunto de instrucciones escritas en un archivo de texto, de forma que estás puedan almacenarse y ejecutarse con posterioridad. Se suelen almacenar con la extensión .R.
Un editor muy recomendado para la creación, ejecución y depuración de estos scripts es RStudio.
? "mean"?? "anova"En RStudio, se puede seleccionar y pulsar F1 para buscar en la ayuda
TRUE o FALSEa <- 2
a## [1] 2
class(a)## [1] "numeric"
b <- 13.6788956789
b## [1] 13.6789
class(b)## [1] "numeric"
print(b, digits = 10)## [1] 13.67889568
Operadores más comunes para datos de tipo numeric
| Suma | + |
| Resta | - |
| Producto | * |
| División | / |
| Potencia | ^ o ** |
| Módulo (residuo) | %% |
| División entera | %/% |
| Comparaciones | == > < >= <= != |
a+b## [1] 15.6789
a-b## [1] -11.6789
a*b## [1] 27.35779
a^b## [1] 13114.69
a/b## [1] 0.1462106
Permite almacenar números enteros. Es el tipo de datos usado para contadores e índices.
n <- as.integer(340000)
class(n)## [1] "integer"
n <- 2L
class(n)## [1] "integer"
Permite almacenar cadenas de texto de longitud indeterminada.
Los literales se indican entre comillas dobles (o simples).
a <- "aaa"
b <- "bbb"
paste(a, b, "hola", sep = ", ") ## [1] "aaa, bbb, hola"
Operadores más comunes para datos de tipo character
| Comparaciones | == > < >= <= != |
Almacena valores que son verdadero (TRUE) o falso (FALSE).
Es el resultado de una comparación.
3 == 2## [1] FALSE
b <- 3 != 2
b## [1] TRUE
Operadores más comunes para datos de tipo logical
| Intersección lógica | & |
| Unión lógica | | |
| Negación lógica | ! |
| Comparaciones | == != |
Por ejemplo…
a <- TRUE
b <- F
a & b # Operador AND## [1] FALSE
a | b # Operador OR## [1] TRUE
!b # Operador NOT## [1] TRUE
Las funciones de conversión de tipos en R se indican como as. seguido del tipo de dato de destino.
as.character(2)## [1] "2"
as.numeric(TRUE)## [1] 1
Para saber si una información o variable es de un tipo determinado puede usarse is. seguido del tipo de dato a comprobar.
R incorpora algunas constantes preestablecidas como
pi## [1] 3.141593
Inf
NaN
NA
NULLTambién existen funciones para evaluar si un dato corresponde a una de estas constantes.
is.na(3)## [1] FALSE
is.null(NULL)## [1] TRUE
is.infinite(-Inf)## [1] TRUE
a <- c(2, 3, 4)
str(a)## num [1:3] 2 3 4
a[2]## [1] 3
Las operaciones se aplican a vectores y devuelven vectores
a + a ## [1] 4 6 8
a > 2.5## [1] FALSE TRUE TRUE
Atención al reciclaje de datos…
a <- c(2, 3, 4)
b <- c(10, 20)
a * b## Warning in a * b: longer object length is not a multiple of shorter object
## length
## [1] 20 60 40
Las funciones en R suelen ser vectoriales
a <- c(2, 3, 4)
abs(sin(a))## [1] 0.9092974 0.1411200 0.7568025
exp(a)## [1] 7.389056 20.085537 54.598150
length(a)## [1] 3
sum(a)## [1] 9
mean(a)## [1] 3
También sd, max, min…
Puede determinarse si un valor está presente en un vector con el operador %in%.
a <- c(2, 3, 4)
3 %in% a## [1] TRUE
c(2,5,3,4) %in% a## [1] TRUE FALSE TRUE TRUE
Existen estructuras especificas para crear vectores secuencia.
1:10## [1] 1 2 3 4 5 6 7 8 9 10
seq(1, 10, by = 2)## [1] 1 3 5 7 9
seq(1, 3, length.out = 5) ## [1] 1.0 1.5 2.0 2.5 3.0
Los vectores se ordenan usando las funciones sort y order.
a <- c(8, 2, 5, 3)
sort(a)## [1] 2 3 5 8
a[order(a)]## [1] 2 3 5 8
a[order(a,decreasing = T)]## [1] 8 5 3 2
Un factor es un vector de cadenas indexado. Normalmente se crea a partir del vector de cadenas de texto.
a <- c("hola", "adeu","hola", "adeu", "adeu", "bye")
b <- as.factor(a)
as.character(b)## [1] "hola" "adeu" "hola" "adeu" "adeu" "bye"
as.numeric(b)## [1] 3 1 3 1 1 2
str(b)## Factor w/ 3 levels "adeu","bye","hola": 3 1 3 1 1 2
La función factor permite la codificación o recodificación manual de un vector.
a <- factor(c(3, 1, 3, 1, 1, 2), labels = c("adeu", "bye", "hola"))
a## [1] hola adeu hola adeu adeu bye
## Levels: adeu bye hola
levels(a)## [1] "adeu" "bye" "hola"
dfA <- data.frame(int = 1:10,
let = sample(letters, 10, replace = TRUE),
ran = rnorm(10))
dfA## int let ran
## 1 1 w 0.3434244
## 2 2 q 0.2019293
## 3 3 i 1.6396732
## 4 4 f 1.8508826
## 5 5 h 0.3068050
## 6 6 n 0.2221603
## 7 7 o 1.8145489
## 8 8 f -0.5816948
## 9 9 e -1.0427837
## 10 10 o -0.7087821
dim(dfA) # Dimensión## [1] 10 3
nrow(dfA) # Número de columnas## [1] 10
ncol(dfA) # Número de filas## [1] 3
str(dfA)## 'data.frame': 10 obs. of 3 variables:
## $ int: int 1 2 3 4 5 6 7 8 9 10
## $ let: Factor w/ 8 levels "e","f","h","i",..: 8 7 4 2 3 5 6 2 1 6
## $ ran: num 0.343 0.202 1.64 1.851 0.307 ...
head(dfA, 3) # Primeros datos, 6 por defecto## int let ran
## 1 1 w 0.3434244
## 2 2 q 0.2019293
## 3 3 i 1.6396732
tail(dfA, 2) # Últimos valores## int let ran
## 9 9 e -1.0427837
## 10 10 o -0.7087821
Para acceder a los datos almacenados en el data frame se usan índices. Las variables de data frame también pueden extraerse usando su nombre.
dfA[2,3]## [1] 0.2019293
dfA[,1]## [1] 1 2 3 4 5 6 7 8 9 10
dfA$let## [1] w q i f h n o f e o
## Levels: e f h i n o q w
Veremos más formas de acceder y usar los datos de una tabla de datos más adelante.
a <- list(2, "2", FALSE)
b <- list(3, "hola", c(2, 3, 4))a## [[1]]
## [1] 2
##
## [[2]]
## [1] "2"
##
## [[3]]
## [1] FALSE
b## [[1]]
## [1] 3
##
## [[2]]
## [1] "hola"
##
## [[3]]
## [1] 2 3 4
length(a)## [1] 3
a[[3]]## [1] FALSE
b[[3]][1]## [1] 2
str(b)## List of 3
## $ : num 3
## $ : chr "hola"
## $ : num [1:3] 2 3 4
notes <- c("Aprovat", "Insuficient", "Notable", "Insuficient",
"Notable", "Excel·lent", "Aprovat")
notes <- factor(notes,
levels = c("Insuficient", "Aprovat",
"Notable", "Excel·lent"),
ordered = TRUE)
str(notes)## Ord.factor w/ 4 levels "Insuficient"<..: 2 1 3 1 3 4 2
levels(notes)## [1] "Insuficient" "Aprovat" "Notable" "Excel·lent"
a <- matrix(c(2, 4, -3, 5), ncol = 2)
a## [,1] [,2]
## [1,] 2 -3
## [2,] 4 5
a[2,2]## [1] 5
a * a # Producto posición por posición## [,1] [,2]
## [1,] 4 9
## [2,] 16 25
a %*% a # Producto matricial## [,1] [,2]
## [1,] -8 -21
## [2,] 28 13
t(a) # Transposición## [,1] [,2]
## [1,] 2 4
## [2,] -3 5
R tiene muchos paquetes para resolver problemas específicos. Para usar un paquete adicional este debe instalarse y cargarse en memoria.
installed.packages()[,1] # Lista los paquetes instalados
(.packages()) # Lista los paquetes en memoriaDos recursos importantes para buscar y identificar paquetes relevantes son https://www.rdocumentation.org/ https://cran.rstudio.com/web/views/
Para instalar un paquete, por ejemplo “nycflights13”
install.packages("nycflights13")Para cargar un paquete en memoria, se usa
library("nycflights13")En RStudio la gestión de paquetes también puede hacerse desde la interfaz del programa
En un script y para garantizar la disponibilidad de un paquete
if(!require("nycflights13")) {
install.packages("nycflights13")
library("nycflights13")
}Para acceder a la documentación de un paquete
help(package="nycflights13")Algunos paquetes tienen información adicional a la que se puede acceder con las funciones vignette() y demo().
La principales funciones de los gráficos en el análisis de datos son
Veremos dos paradigmas distintos para la creación de gráficos en R.
base.ggplot2base RLos gráficos en base R se construyen a partir de datos en vectores y usando funciones específicas en función del gráfico a realizar.
Para ejemplificar algunas de estas funciones, usaremos el conjunto de datos mtcars, que forma de los paquetes básicos de R.
str(mtcars)## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...
head(mtcars, 10)## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
Veremos como crear
base R – gráfico de dispersiónplot(mtcars$disp,mtcars$hp)base R – histogramahist(mtcars$disp)Añadiendo una curva de densidad…
hist(mtcars$disp,breaks=20,freq=FALSE)
lines(density(mtcars$disp))base R – diagrama de barrasbarplot(table(as.factor(mtcars$cyl)))base R – diagrama de cajaboxplot(mtcars$disp)
points(mean(mtcars$disp))La gramática de gráficos es una aproximación teórica al estudio de los componentes de un gráfico. De acuerdo con este análisis, un gráfico se puede construir mediante la especificación de un conjunto de capas y componentes que definen los datos, la asociación de los mismos a aspectos del gráfico, la especificación de la relación entre los valores de las variables de los datos con las del gráfico, la estructura geométrica del gráfico…
Wilkinson, L. (2006). The grammar of graphics. Springer Science & Business Media.
ggplot2ggplot2 es un paquete de R que implementa de la gramática de gráficos.
Wickham, H. (2010). A layered grammar of graphics. Journal of Computational and Graphical Statistics, 19(1), 3-28.
Referencias:
ggplot2 forma parte del paquete tidyverse (aunque también puede instalarse y cargarse autónomamente).
if(!require("tidyverse")) {
install.packages("tidyverse")
library("tidyverse")
}ggplot2 – capas y elementos del gráficoEn ggplot2 cada elemento gráfico que representa un conjunto de datos constituye una capa. Una o más capas constituyen un gráfico.
Cada capa queda definida mediante la especificación de sus elementos. Los principales son
En ggplot2, los gráficos constituyen un objeto de R y se construyen de forma aditiva.
Por ejemplo,
grafico <- ggplot(data = anscombe,
mapping = aes(x = x1, y = y1)) # Datos y mapeado estético
grafico <- grafico + geom_point() # Geometría
graficoggplot2 – datosEn ggplot2, el elemento data (datos) se introduce como primer argumento de la función ggplot. Debe corresponder a una tabla de datos o un tipo de datos convertible a tabla de datos.
grafico <- ggplot(data = anscombe,ggplot2 – mapeado estéticoEl mapping (mapeado estético) corresponde al establecimiento de relaciones entre variables de los datos y variables del gráfico. Es el segundo argumento de la función ggploty debe crearse con al función de apoyo aes.
mapping = aes(x = x1, y = y1))Para variables cuantitativas, los mapeados más comunes corresponden a
x, y…sizecolor, fillPara variable cualitativas, los mapeados más frecuentes son
x, y…color, fillshapeggplot2 – geometríasLas geometrías (geom_) indican la forma que debe tener el gráfico, es decir, cómo se articulan las variables del gráfico. Se añaden al gráfico sumándose al objeto creado por ggplot.
grafico <- grafico + geom_point()Son geometrías de uso común
geom_pointgeom_line, geom_vline, geom_hlinegeom_bargeom_histogramgeom_boxplotUn resumen de las geometrías y su relación con las variables del gráfico que reconoce cada una de ellas figura en https://github.com/rstudio/cheatsheets/raw/master/data-visualization-2.1.pdf.
ggplot2Veremos cómo crear en ggplot2 los gráficos más habituales, añadiendo algunas consideraciones para aquellos casos donde los gráficos realizados con base tienen prestaciones insuficientes.
De forma general, los gráficos en ggplot2 se construyen a partir de tablas de datos (data frames), de los cuales se seleccionan las variables a representar.
Usaremos 1000 datos del conjunto de datos diamonds para crear los distintos ejemplos.
diaM <- diamonds[sample(1:nrow(diamonds),1000),]
str(diaM)## Classes 'tbl_df', 'tbl' and 'data.frame': 1000 obs. of 10 variables:
## $ carat : num 0.41 1.16 0.3 0.41 0.24 1.11 0.51 1.31 0.41 1.16 ...
## $ cut : Ord.factor w/ 5 levels "Fair"<"Good"<..: 2 4 4 5 3 5 5 5 5 5 ...
## $ color : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 2 4 1 5 2 4 4 4 2 6 ...
## $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 3 4 6 8 6 8 6 4 5 2 ...
## $ depth : num 63.5 59.9 62.4 62.9 61.6 61.2 61.8 61.9 62.2 62.4 ...
## $ table : num 57 60 58 54 59 54 56 56 55 55 ...
## $ price : int 755 7203 965 1243 552 10053 1875 9085 1030 3800 ...
## $ x : num 4.69 6.88 4.3 4.74 3.98 6.71 5.1 7.08 4.79 6.76 ...
## $ y : num 4.73 6.85 4.26 4.71 4.01 6.73 5.16 7.05 4.76 6.69 ...
## $ z : num 2.99 4.11 2.67 2.97 2.46 4.11 3.17 4.37 2.97 4.2 ...
ggplot2 – gráfico de dispersiónggplot(diaM, aes(x=carat,y=price)) + geom_point()Añadiendo una tercera variable (cut) y modificando algunos aspectos de formato…
ggplot(diaM, aes(x=carat,y=price,color=cut)) +
geom_point(alpha=.8,shape=21,size=3)Añadiendo líneas de tendencia…
ggplot(diaM, aes(x=carat,y=price,color=cut)) +
geom_point(alpha=.8,shape=21,size=3) +
geom_smooth(method="lm",se=FALSE)ggplot2 – histogramaggplot(diaM, aes(x=price)) + geom_histogram(binwidth=1000)Y en función del corte…
ggplot(diaM, aes(x=price,fill=cut)) +
geom_histogram(position='dodge',binwidth=1000)En frecuencias relativas (por grupo)…
ggplot(diaM, aes(x=price,y=..density..,fill=cut)) +
geom_histogram(position='dodge',binwidth=1000)Quizás funcione mejor un gráfico de densidades…
ggplot(diaM, aes(x=price,fill=cut)) +
geom_density(alpha=.3)ggplot2 – diagrama de barrasggplot(diaM, aes(x=clarity)) + geom_bar()En función de la claridad…
ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar()Para comparar entre frecuencias absolutas, funcionan mejor las barras separadas.
ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar(position="dodge")Para comparar entre frecuencias relativas acumuladas, son mejores las barras apiladas en frecuencia relativa (para cada clase).
ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar(position="fill")Los diagramas de barras también pueden crearse a partir de tablas de datos agrupados. En este caso, debe indicarse qué variable es la y e incluir stat="identity" en el geom_bar.
ggplot2 – diagrama de cajaggplot(diaM, aes(x=1, y=price)) + geom_boxplot()NOTA: Para el geom_boxplot se requieren dos variables. Si no hay variable independiente se puede incluir x=1.
Y en función del corte…
ggplot(diaM, aes(x=cut, y=price)) + geom_boxplot()El gráfico se puede mejorar mostrando todos los puntos, con una posición aleatorizada.
ggplot(diaM, aes(x=cut, y=price)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(shape = 21, alpha=.5,height=0,width=.2)O incluyendo un violin plot y un punto para la media…
medias <- diaM %>% group_by(cut) %>%
summarise(price=mean(price))
ggplot(diaM, aes(x=cut, y=price)) +
geom_violin() +
geom_boxplot(outlier.shape = NA, width = 0.1) +
geom_point(data=medias,shape=3)ggplot2 – elementos adicionalesAdemás de los elementos ya presentados (datos, mapeado estético y geometrías), otros elementos de ggplot2 permiten controlar aspectos adicionales del gráfico, por ejemplo
scale_...: controlan los aspectos relativos a la presentación de las variables del gráfico,coord_...: establecen el sistema de coordenadas usado en la geometría,labs: establece los títulos del gráfico,theme, theme_...: controlan la part del gráfico que no corresponde a datos (non-data ink), yfacet_: permite la creació de secuencia de gráficos en función de una o dos variablesUn ejemplo para terminar…
ggplot(diaM, aes(x=carat, y = price, shape = cut, col = clarity)) +
geom_point(alpha=.6) +
scale_x_continuous(breaks=1:3) +
scale_y_continuous(trans="log10") +
scale_color_brewer(palette="Spectral")+
facet_grid(cut ~ clarity) +
theme_bw() +
theme(legend.position = "none",text = element_text(size=10))La estadística, como parte de la matemática que se ocupa de la recolección, análisis e interpretación de datos, tiene dos funciones principales
En este repaso de la estadística básica, usaremos (de nuevo) el conjunto de datos mtcars.
str(mtcars)## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...
head(mtcars, 10)## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
## Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
## Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
A menudo, el primer interés que tendremos es el de describir un conjunto de datos correspondientes a una sola variable. Esta es la situación con la que empezaremos.
¿Cómo se describe un conjunto de datos univariante?
Normalmente se consideran los siguientes aspectos
Tamaño de la muestra o número de datos
length(mtcars$disp)## [1] 32
n <- sum(!is.na(mtcars$disp))
n## [1] 32
Frecuencias absolutas y relativas
Para una variable cuantitativa (numeric)…
frec_abs <- table(cut(mtcars$disp,
breaks=c(0,100,200,300,400,500,Inf)))
frec_abs##
## (0,100] (100,200] (200,300] (300,400] (400,500] (500,Inf]
## 5 11 5 8 3 0
frec_rel <- frec_abs / n
frec_rel##
## (0,100] (100,200] (200,300] (300,400] (400,500] (500,Inf]
## 0.15625 0.34375 0.15625 0.25000 0.09375 0.00000
Cuantiles
quantile(mtcars$disp,.1)## 10%
## 80.61
quantile(mtcars$disp,1:4*.2)## 20% 40% 60% 80%
## 120.14 160.00 275.80 350.80
mean(mtcars$disp)## [1] 230.7219
median(mtcars$disp)## [1] 196.3
var(mtcars$disp) # el denominador es n-1## [1] 15360.8
sd(mtcars$disp)## [1] 123.9387
range(mtcars$disp)## [1] 71.1 472.0
diff(range(mtcars$disp))## [1] 400.9
IQR(mtcars$disp)## [1] 205.175
mad(mtcars$disp) # mean absolute deviation## [1] 140.4764
lmin <- quantile(mtcars$disp,.25) - 1.5*IQR(mtcars$disp)
lmax <- quantile(mtcars$disp,.75) + 1.5*IQR(mtcars$disp)
mtcars$disp[mtcars$disp > lmax | mtcars$disp < lmin]## numeric(0)
lmin <- quantile(mtcars$qsec,.25) - 1.5*IQR(mtcars$qsec)
lmax <- quantile(mtcars$qsec,.75) + 1.5*IQR(mtcars$qsec)
mtcars$qsec[mtcars$qsec > lmax | mtcars$qsec < lmin]## [1] 22.9
Gráficamente, son útiles para describir un conjunto de datos
boxplot(mtcars$qsec)
hist(mtcars$qsec)barplot(table(factor(mtcars$cyl)))Para el estudio y descripción de la relación entre dos variables, las técnicas más habituales incluyen
Es el concepto equivalente a las tablas de frecuencias en el caso univariante.
Para variables cualitativas…
table(factor(mtcars$am),factor(mtcars$cyl))##
## 4 6 8
## 0 3 4 12
## 1 8 3 2
Para variables cuantitativas…
table(cut(mtcars$disp,breaks = seq(0,500,by=100)),
cut(mtcars$qsec,breaks = seq(13,25,by=3)))##
## (13,16] (16,19] (19,22] (22,25]
## (0,100] 0 3 2 0
## (100,200] 1 7 2 1
## (200,300] 0 3 2 0
## (300,400] 4 4 0 0
## (400,500] 0 3 0 0
Producto-momento de Pearson
cor(mtcars$disp,mtcars$qsec)## [1] -0.4336979
Coeficiente de correlación de Spearman
cor(mtcars$disp,mtcars$qsec, method = "spearman")## [1] -0.4597818
plot(mtcars$disp,mtcars$qsec)La distribución de probabilidad de un número aleatorio corresponde a la abstracción teórica de la densidad de los datos de una conjunto. Representa la densidad que se obtendría cunado se generasen infinitos números aleatorios de acuerdos a un mismo procedimiento o experimento.
Un gran número de distribuciones experimentales tienen modelos teóricos con los que se relacionan. Los tres más comunes son la distribución normal (para datos de variables cuantitativas continuas), la distribución uniforme (para cualquier tipo de datos) y la distribución binomial (para el número de veces que se produce un suceso de una determinada probabilidad).
En R, todas las distribuiciones comparten el mismo sistema de funciones
r<dist>: genera aleatorios deacuerdo con una distribución de probabilidadd<dist>: devuelve la densidad de probablidad para un valor de variableq<dist>: devuelve el valor de la variable tal que la probablidadp<dist>: devuelve la probablidad que un valor sea inferior o igual a xPor ejemplo, para la distribución normal estándar…
rnorm(10)## [1] -0.5040593 -0.4025292 -1.2126037 -0.7030878 -0.2707858 0.2570524
## [7] 0.6062162 -0.1483702 -1.3017443 0.4667901
dnorm(0)## [1] 0.3989423
qnorm(.95)## [1] 1.644854
pnorm(1.64)## [1] 0.9494974
df <- data.frame(x = rnorm(1000, mean = 3, sd = 1))
dfT <-data.frame(x = seq(0,6,length.out=101),
y = dnorm(seq(0,6,length.out=101),mean=3,sd=1))
ggplot(df, aes(x = x, y=..density..)) +
geom_histogram(binwidth = .1,fill="#dddddd",col="black") +
geom_density() +
geom_line(data=dfT,aes(x=x,y=y),col="red")+
theme_bw()pnorm(5,mean = 3,sd = 1)## [1] 0.9772499
qnorm(.98,mean = 3,sd = 1)## [1] 5.053749
df <- data.frame(x = runif(1000, min = 10, max = 20))
dfT <-data.frame(x = seq(10,20,length.out=101),
y = dunif(seq(10,20,length.out=101), min=10, max=20))
ggplot(df, aes(x = x, y=..density..)) +
geom_histogram(binwidth = .5,fill="#dddddd",col="black") +
geom_density() +
geom_line(data=dfT,aes(x=x,y=y),col="red")+
theme_bw()punif(12, min=10, max=20)## [1] 0.2
qunif(.90, min=10, max=20)## [1] 19
df <- data.frame(x = rbinom(100,5,prob=0.5))
df <- df %>% group_by(x) %>%
summarise(y = n()/nrow(df))
dfT <- data.frame(x = 0:5,
y = dbinom(0:5,size=5,prob=0.5))
dfJ <- rbind(cbind(teo=FALSE,df),
cbind(teo=TRUE,dfT))
ggplot(dfJ,aes(x=x,y=y,fill=teo)) +
geom_bar(position="dodge", stat="identity",
col="black") +
scale_fill_manual(values=c("#dddddd","red")) +
theme_bw() +
theme(legend.position = "none")pbinom(2,5,prob=0.5)## [1] 0.5
qbinom(.5,5,.5)## [1] 2
R incorpora otras muchas distribuciones de probabilidad que pueden consultarse en la página correspondiente de la ayuda.
? "Distributions"En términos estadísticos, la inferencia consiste en la obtención de información sobre la población –el conjunto de los valores existentes– a partir de una muestra representativa de la misma.
Si los datos que tenemos, no constituyen una muestra representativa cualquier información que queramos extrapolar a la población estará sesgada, es decir, desplazada de su valor real.
Los resultados de la inferencia se suelen mostrar mediante una de las dos formas siguientes:
Presentaremos a continuación, las técnicas más comúnmente usadas para
Para contrastar la hipótesis que un conjunto de datos puede representarse mediante un determinada distribución teórica, las pruebas más comunes son la prueba de bondad de ajuste de chi cuadrado (chisq.test) y la prueba de Kolmogorov-Smirnov (ks.test). En ambos casos, los argumentos son las frecuencias absolutas de cada clase para la distribución experimental y las probabilidades esperadas de acuerdo con la distribución teórica o de referencia.
Ejemplo
Se ha tirado cien veces un dado de 10 caras, y se han obtenido cada una de las caras las veces que se muestran en la tabla siguiente.
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| 9 | 6 | 7 | 15 | 11 | 11 | 9 | 13 | 9 | 10 |
Si el dado fuese normal, todas las caras tendrían la misma probabilidad de aparecer. ¿Lo es?
obs <- c(9, 6, 7, 15, 11, 11, 9, 13, 9, 10)
exp <- rep(.1, 10)
chisq.test(x = obs, p = exp)##
## Chi-squared test for given probabilities
##
## data: obs
## X-squared = 6.4, df = 9, p-value = 0.6993
Gráficamente, y aunque sin la misma objetividad, la comparación de distribuciones también suele hacerse mediante el gráfico de cuantiles –en R, es la función qqplot.
Si ambos conjuntos de datos corresponden a la misma distribución, sus cuantiles deben quedar alineados junto a la diagonal del gráfico (especialemente en la región central del mismo).
exp <- rep(1:10, obs)
teo <- rep(1:10, 10)
qqplot(x = exp, y = quantile(teo, (1:length(exp))/length(exp)))
qqline(y = exp, distribution = function(x) quantile(teo,x),
probs= c(1/length(exp),1))Para comprobar si un conjunto de datos puede proceder de una población aleatoria de distribución normal, existen pruebas estadísticas más específicas. Entre las muchas pruebas existentes, son de uso común las pruebas de Shapiro-Wilk y Anderson-Darling.
x1 <- rnorm(20, mean = 4, sd = 5)
x2 <- rbeta(20, shape1 = 5, shape2 = .5, ncp = 4)
shapiro.test(x1)##
## Shapiro-Wilk normality test
##
## data: x1
## W = 0.97713, p-value = 0.8919
shapiro.test(x2)##
## Shapiro-Wilk normality test
##
## data: x2
## W = 0.59299, p-value = 2.493e-06
if(!require("nortest")) {
install.packages("nortest")
library("nortest")
}ad.test(x1)##
## Anderson-Darling normality test
##
## data: x1
## A = 0.24361, p-value = 0.7311
ad.test(x2)##
## Anderson-Darling normality test
##
## data: x2
## A = 3.4539, p-value = 5.417e-09
Gráficamente se puede usar el gráfico de cuantiles comentado anteriormente.
qqnorm(x1)
qqline(x1)
qqnorm(x2)
qqline(x2)Las inferencias más habituales respecto a la dispersión de una población constituyen los siguientes casos
Cuando los datos estan normalmente distribuidos, el intervalo de confianza de la varianza se calcula a partir de la distribución de chi cuadrado.
x <- rnorm(50, mean = 0, sd = 2)
df <- length(x) - 1
lower <- var(x) * df / qchisq(1 - 0.05/2, df)
upper <- var(x) * df / qchisq(0.05/2, df)
c(lower = lower, var = var(x), upper = upper)## lower var upper
## 3.027845 4.339236 6.738175
Para obtener el intervalo de confianza para la desviación estándar, cuando los datos estan normalmente distribuidos, se calcula la desviación estándar a partir de las varainzas calculadas más arriba.
c(lower = sqrt(lower), sd = sd(x), upper = sqrt(upper))## lower sd upper
## 1.740070 2.083083 2.595800
Para la comparación de la dispersión de las poblaciones para dos conjuntos de datos, cuando los datos estan normalmente distribuidos, se lleva a cabo con un prueba de F –función var.test en R–.
x <- rnorm(50, mean = 0, sd = 2)
y <- rnorm(30, mean = 1, sd = 1)
var.test(x, y)##
## F test to compare two variances
##
## data: x and y
## F = 4.4128, num df = 49, denom df = 29, p-value = 5.781e-05
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 2.217108 8.302372
## sample estimates:
## ratio of variances
## 4.41283
Si los datos no estan normalmente distribuidos, la comparación de dispersiones puede hacerse mediante el test de Ansari.
x <- rlnorm(50, meanlog = 2, sdlog = 1)
y <- rlnorm(30, meanlog = 2, sdlog = .2)
ansari.test(x, y)##
## Ansari-Bradley test
##
## data: x and y
## AB = 823, p-value = 5.921e-05
## alternative hypothesis: true ratio of scales is not equal to 1
Para comparar las dispersiones de más de un conjunto de datos para los cuáles se pueda asumir normalidad, se usa la prueba de Bartlett –bartlett.test en R–.
x1 <- round(rnorm(20, mean = 1, sd = 2),1)
x2 <- round(rnorm(20, mean = 3, sd = 2),1)
x3 <- round(rnorm(20, mean = 5, sd = 2),1)
list(x1,x2,x3)
bartlett.test(list(x1,x2,x3))## [[1]]
## [1] 4.6 3.6 1.1 0.1 -2.8 1.3 3.9 1.7 -2.6 4.0 0.1 4.8 1.5 0.7
## [15] -2.9 1.3 3.1 -3.1 1.4 -0.5
##
## [[2]]
## [1] 1.3 3.0 2.8 2.8 -2.0 0.5 4.4 4.2 3.6 7.2 4.2 0.2 0.7 2.9
## [15] 6.5 5.2 4.2 3.6 2.0 4.3
##
## [[3]]
## [1] 4.7 3.2 7.7 4.3 6.9 6.0 1.3 4.4 4.0 5.6 4.3 0.9 4.0 8.5 2.2 3.3 5.4
## [18] 6.8 2.9 4.0
##
## Bartlett test of homogeneity of variances
##
## data: list(x1, x2, x3)
## Bartlett's K-squared = 0.97116, df = 2, p-value = 0.6153
Si los datos no cumplen con la hipótesis de normalidad, en lugar de la prueba de Bartlett se usa la prueba de Levene o la de Brown–Forsythe –leveneTest en el paquete car de R–.
Gráficamente y aunque de forma menos objetiva, la comparación de dispersiones puede hacerse usando gráficos de caja, bien directamente, bien centrando los datos.
x1 <- round(rnorm(20, mean = 1, sd = 2),1)
x2 <- round(rnorm(20, mean = 3, sd = 2),1)
x3 <- round(rnorm(20, mean = 5, sd = 2),1)
xs <- data.frame(
group = factor(c(rep(1, length(x1)),rep(2, length(x2)),
rep(3, length(x3)))),
y = c(x1,x2,x3),
centered = c(x1 - median(x1),x2 - median(x2),x3 - median(x3)))
ggplot(xs, aes(x=group,y=y)) + geom_boxplot() +
coord_flip() + theme_bw()
ggplot(xs, aes(x=group,y=centered)) + geom_boxplot() +
coord_flip() + theme_bw()Las principales inferencias respecto a la localización –o tendencia central– de los datos corresponden a
Para establecer el intervalo de confianza de la media, cuando los datos están normalmente distribuidos, se usa la distribución t. El intervalo puede calcularse a partir de la función qt o visualizarlo como resultado de la función t.test.
x <- rnorm(10, mean = 2, sd = .5)
x## [1] 1.975816 1.157617 2.268955 3.450344 2.088544 1.991210 1.575897
## [8] 1.727211 2.054207 2.134882
t.test(x)##
## One Sample t-test
##
## data: x
## t = 10.915, df = 9, p-value = 1.718e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 1.619169 2.465768
## sample estimates:
## mean of x
## 2.042468
Si los datos no están normalmente distribuidos, la función wilcox.test ofrece entre sus resultados el intervalo de confianza para la mediana.
x <- rnorm(10, mean = 3, sd = .5) ^ 3
x## [1] 12.84760 44.27621 34.72355 32.56462 35.15335 41.52064 42.10572
## [8] 25.93625 38.23794 61.32878
wilcox.test(x, conf.int = TRUE)##
## Wilcoxon signed rank test
##
## data: x
## V = 55, p-value = 0.001953
## alternative hypothesis: true location is not equal to 0
## 95 percent confidence interval:
## 27.47666 44.27621
## sample estimates:
## (pseudo)median
## 37.08819
Las dos funciones usadas para obtener los intervalos de confianza permiten comparar la tendencia central de un conjunto de datos con un valor preestablecido. Como se ha indicado antes, la prueba de t require normalidad de los datos; la prueba de Wilcoxon, no.
x <- rnorm(10, mean = 2, sd = .5)
x## [1] 1.273417 1.719438 2.449046 2.087944 2.312369 2.219420 1.094407
## [8] 1.654159 2.895084 2.289877
t.test(x, mu = 1.5)##
## One Sample t-test
##
## data: x
## t = 2.8388, df = 9, p-value = 0.01944
## alternative hypothesis: true mean is not equal to 1.5
## 95 percent confidence interval:
## 1.601462 2.397570
## sample estimates:
## mean of x
## 1.999516
wilcox.test(x ^ 3, mu = 8)##
## Wilcoxon signed rank test
##
## data: x^3
## V = 33, p-value = 0.625
## alternative hypothesis: true location is not equal to 8
Para comparar la localización de dos conjuntos de datos, deben tenerse en cuenta los siguientes aspectos
Si los datos están asociados, es decir que existen pares de datos en ambas muestras tales que comparten fuentes de variación, entonces
t.test o bien usar directamente la función t.test con la opción paired = TRUE,wilcox.test sobre la resta o con la opción paired = TRUE en R–.x <- rnorm(10, mean = 3, sd = 1)
y <- rnorm(10, mean = 3.5, sd = 1)
list(x = x,y = y)## $x
## [1] 3.345022 4.567567 2.504564 2.523306 5.007635 2.739927 3.147583
## [8] 4.118539 2.384038 3.360721
##
## $y
## [1] 1.891288 4.138754 2.288343 5.170342 4.556190 4.928725 5.075908
## [8] 2.936908 3.863685 6.112612
t.test(x - y)##
## One Sample t-test
##
## data: x - y
## t = -1.4082, df = 9, p-value = 0.1927
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.8932433 0.4404728
## sample estimates:
## mean of x
## -0.7263852
t.test(x, y, paired = TRUE)##
## Paired t-test
##
## data: x and y
## t = -1.4082, df = 9, p-value = 0.1927
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.8932433 0.4404728
## sample estimates:
## mean of the differences
## -0.7263852
x <- rlnorm(10, meanlog = 3, sdlog = 1)
y <- rlnorm(10, meanlog = 3.5, sdlog = 1)
list(x = x,y = y)## $x
## [1] 16.275599 118.552003 3.979148 4.856603 92.009866 10.225033
## [7] 18.204673 115.666602 9.336026 18.787804
##
## $y
## [1] 22.68449 21.23973 17.02835 99.80393 16.41797 43.08411 16.15939
## [8] 25.74315 85.43956 165.29406
wilcox.test(x - y)##
## Wilcoxon signed rank test
##
## data: x - y
## V = 22, p-value = 0.625
## alternative hypothesis: true location is not equal to 0
wilcox.test(x, y, paired = TRUE)##
## Wilcoxon signed rank test
##
## data: x and y
## V = 22, p-value = 0.625
## alternative hypothesis: true location shift is not equal to 0
Si los datos no están asociados, entonces
t.test; la opción var.equal permite escoger entre la prueba calculada con la varianza de los datos, y la aproximación de Welch, que se usa cuando la varianzas no son iguales.wilcox.test en R–.x <- rnorm(10, mean = 3, sd = 1)
y <- rnorm(14, mean = 3.5, sd = 2)
list(x = x,y = y)## $x
## [1] 2.756588 4.963940 2.125890 2.265344 2.790045 2.927530 2.223230
## [8] 2.508153 3.366607 1.868743
##
## $y
## [1] -0.03273222 0.37884725 2.53149830 4.33183571 3.04812589
## [6] 2.78806269 2.21554483 0.85147586 6.08607800 2.70657752
## [11] 3.09383981 2.51752542 2.33729053 4.52201226
t.test(x, y, var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: x and y
## t = 0.21253, df = 20.832, p-value = 0.8338
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.9659723 1.1857600
## sample estimates:
## mean of x mean of y
## 2.779607 2.669713
x <- rlnorm(10, meanlog = 3, sdlog = 1)
y <- rlnorm(14, meanlog = 3.5, sdlog = 2)
list(x = x,y = y)## $x
## [1] 16.420154 6.641600 19.814867 4.027625 25.761614 11.242345 16.520460
## [8] 51.367157 29.315267 19.482693
##
## $y
## [1] 442.209234 185.794834 325.379738 31.431309 99.299656
## [6] 543.129873 1125.963140 66.345844 33.067248 214.198444
## [11] 54.738450 1.569378 957.049102 99.277018
wilcox.test(x, y)##
## Wilcoxon rank sum test
##
## data: x and y
## W = 12, p-value = 0.0002743
## alternative hypothesis: true location shift is not equal to 0
Para la comparación de tres o más grupos, las técnicas más habituales son
aov en R–;krukal.test–.x <- c(rnorm(10, mean = 3, sd = 1),
rnorm(8, mean = 3.2, sd = 1),
rnorm(10, mean = 4, sd = 1))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
test <- aov(x ~ g)
summary(test)## Df Sum Sq Mean Sq F value Pr(>F)
## g 2 5.094 2.5469 2.824 0.0784 .
## Residuals 25 22.546 0.9019
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
x <- c(runif(10, min = 2, max = 4),
runif(8, min = 1.5, max = 3.5),
runif(10, min = 3, max = 5))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
kruskal.test(x ~ g)##
## Kruskal-Wallis rank sum test
##
## data: x by g
## Kruskal-Wallis chi-squared = 16.925, df = 2, p-value = 0.0002113
Tanto el análisis de la varianza como la prueba de Kruskal-Wallis solo permiten contrastar si algún par de grupos presentan valores de localización distintos, pero no permiten identificar cuál o cuáles.
Para identificar qué pares de grupos presentan diferencias estadísticamente significativos se usan pruebas post-hoc
TukeyHSD en R,dunn.test en el paquete dunn.test de R–Gráficamente y aunque con menos objetividad, la comparación de la loclización puede hacerse usando gráficos de caja.
x <- c(rnorm(10, mean = 3, sd = 1),
rnorm(8, mean = 3.2, sd = 1),
rnorm(10, mean = 4, sd = 1))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
xs <- data.frame(y = x, group = g)
ggplot(xs, aes(x=group,y=y)) + geom_boxplot() +
coord_flip() + theme_bw()El ajuste de modelos se expresa de forma habitual en R mediante un objeto formula.
En esta notación, se indica la relación entre variables mediante una expresión de tres términos donde la variable dependiente se indica a la izquierda y las variables dependientes a la derecha.
form1 <- y ~ x # recta
class(form1)## [1] "formula"
form1 <- y ~ log(x) # logaritmo
form1 <- y ~ poly(x,4) # polinomio de grado 4
form1 <- y ~ x + 0 # recta que pasa por el origen
form1 <- y ~ I(x^.5) # raíz cuadradaEl ajuste de un modelo lineal por mínimos cuadrados ordinarios se hace en R mediante la función lm.
fit1 <- lm(mpg ~ wt, data=mtcars)
summary(fit1)##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10
Una vez ajustado el modelo, este puede usarse para predecir nuevos datos (o calcular los valores ajustados) mediante la función predict.
df <- data.frame(x=mtcars$wt, y=mtcars$mpg,
predict(fit1, newdata=mtcars, interval="prediction"))
head(df)## x y fit lwr upr
## Mazda RX4 2.620 21.0 23.28261 16.92894 29.63628
## Mazda RX4 Wag 2.875 21.0 21.91977 15.59072 28.24882
## Datsun 710 2.320 22.8 24.88595 18.48644 31.28546
## Hornet 4 Drive 3.215 21.4 20.10265 13.78568 26.41962
## Hornet Sportabout 3.440 18.7 18.90014 12.57806 25.22223
## Valiant 3.460 18.1 18.79325 12.47021 25.11630
El modelo se puede visualizar gráficamente a partir de estas predicciones.
df2 <- data.frame(
wt = seq(min(mtcars$wt), max(mtcars$wt), length.out = 201),
predict(fit1,
newdata = data.frame(wt = seq(min(mtcars$wt), max(mtcars$wt),
length.out = 201)),
interval="prediction"))
ggplot(mtcars, aes(x=wt, y=mpg)) +
geom_ribbon(data=df2, mapping = aes(x=wt, ymin=lwr,ymax=upr),
inherit.aes = FALSE, fill="#dddddd", alpha=0.8)+
geom_point(shape=3) +
geom_line(data=df2, mapping = aes(x=wt, y=fit)) +
theme_bw()Para comprobar la adecuación del modelo y del método de ajuste usado (OLS) deben analizarse los residuales.
Estos deben ser aleatorios, estar normalmente distribuidos, mostrar homocedasticidad y no mostrar evidencia de puntos especialmente influyentes.
Estas comprobaciones pueden hacerse a partir de los datos de los residuales y las pruebas estadísticas oportunas.
fit1$residuals## Mazda RX4 Mazda RX4 Wag Datsun 710
## -2.2826106 -0.9197704 -2.0859521
## Hornet 4 Drive Hornet Sportabout Valiant
## 1.2973499 -0.2001440 -0.6932545
## Duster 360 Merc 240D Merc 230
## -3.9053627 4.1637381 2.3499593
## Merc 280 Merc 280C Merc 450SE
## 0.2998560 -1.1001440 0.8668731
## Merc 450SL Merc 450SLC Cadillac Fleetwood
## -0.0502472 -1.8830236 1.1733496
## Lincoln Continental Chrysler Imperial Fiat 128
## 2.1032876 5.9810744 6.8727113
## Honda Civic Toyota Corolla Toyota Corona
## 1.7461954 6.4219792 -2.6110037
## Dodge Challenger AMC Javelin Camaro Z28
## -2.9725862 -3.7268663 -3.4623553
## Pontiac Firebird Fiat X1-9 Porsche 914-2
## 2.4643670 0.3564263 0.1520430
## Lotus Europa Ford Pantera L Ferrari Dino
## 1.2010593 -4.5431513 -2.7809399
## Maserati Bora Volvo 142E
## -3.2053627 -1.0274952
Gráficamente, los mismas comprobaciones pueden hacerse a partir de la representación gráfica del modelo.
plot(fit1, which=1:6)Estas mismas ideas son la base de la aplicación en R de las técnicas de ajuste multilineal (lm), ajuste lineal generalizado (glm) y ajuste no lineal por mínimos cuadrados (nls).
A menudo es necesario leer y guardar los datos en algún formato tal que permita la importación y exportación de los mismos y su intercambio con otros programas o entornos.
Aunque puede importar y exportar archivos de muy distintos formatos, solo comentaremos como trabajar con
Para los distintos ejemplos usaremos el conjunto de datos flights del paquete nycflights13.
if(!require("nycflights13")) {
install.packages("nycflights13")
library("nycflights13")
}str(flights)## Classes 'tbl_df', 'tbl' and 'data.frame': 336776 obs. of 19 variables:
## $ year : int 2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...
## $ month : int 1 1 1 1 1 1 1 1 1 1 ...
## $ day : int 1 1 1 1 1 1 1 1 1 1 ...
## $ dep_time : int 517 533 542 544 554 554 555 557 557 558 ...
## $ sched_dep_time: int 515 529 540 545 600 558 600 600 600 600 ...
## $ dep_delay : num 2 4 2 -1 -6 -4 -5 -3 -3 -2 ...
## $ arr_time : int 830 850 923 1004 812 740 913 709 838 753 ...
## $ sched_arr_time: int 819 830 850 1022 837 728 854 723 846 745 ...
## $ arr_delay : num 11 20 33 -18 -25 12 19 -14 -8 8 ...
## $ carrier : chr "UA" "UA" "AA" "B6" ...
## $ flight : int 1545 1714 1141 725 461 1696 507 5708 79 301 ...
## $ tailnum : chr "N14228" "N24211" "N619AA" "N804JB" ...
## $ origin : chr "EWR" "LGA" "JFK" "JFK" ...
## $ dest : chr "IAH" "IAH" "MIA" "BQN" ...
## $ air_time : num 227 227 160 183 116 150 158 53 140 138 ...
## $ distance : num 1400 1416 1089 1576 762 ...
## $ hour : num 5 5 5 5 6 5 6 6 6 6 ...
## $ minute : num 15 29 40 45 0 58 0 0 0 0 ...
## $ time_hour : POSIXct, format: "2013-01-01 05:00:00" "2013-01-01 05:00:00" ...
Partiremos de un subconjunto de flights, para ello empezamos segmentando el data.frame.
fl_ny2ws <- flights[flights$dest %in% c("IAD","BWI"),
c("origin","dest","carrier","arr_delay","air_time")]
head(fl_ny2ws)## # A tibble: 6 x 5
## origin dest carrier arr_delay air_time
## <chr> <chr> <chr> <dbl> <dbl>
## 1 LGA IAD EV -14 53
## 2 LGA BWI WN -19 40
## 3 EWR IAD EV 12 52
## 4 JFK BWI MQ 851 41
## 5 JFK IAD B6 14 52
## 6 LGA IAD EV 5 59
Para escribir un archivo de texto (delimitado) desde un data.frame se usa la función write.table o cualquiera de sus derivadas. Para leer, la función es read.table.
write.table(fl_ny2ws, file="fl_ny2ws.csv", sep=",", dec=".",
quote=TRUE, fileEncoding="UTF-8", row.names=FALSE)fl_ny2ws <- read.table("fl_ny2ws.csv", sep=",", dec=".",
quote="\"", fileEncoding="UTF-8", header=TRUE)Si se desea o dispone de un archivo delimitado por tabulaciones, el carácter tabulador se indica como \t.
Para leer y guardar datos en el formato propio de R, se usan las funciones save y load. Estas permiten almacenar y recuperar cualquier conjunto de variables del entorno de trabajo. Al recuperarlas se recuperan con el mismo nombre con el que se almacenaron.
save(fl_ny2ws, file="fl_ny2ws.rda")print(load("fl_ny2ws.rda")) #print muestra el nombre de los objetosCom se ha viso anteriormente, el data frame es el tipo de dato más usado para almacenar tablas de datos.
A menudo, para poder realizar gráficos y/o aplicar distintos procedimientos estadísticos, es necesario manipular la tabla de datos. A esto dedicaremos este apartado.
Entre las operaciones habituales que haremos sobre una tabla de datos, hay
Terminaremos este bloque introduciendo los funciones básicas del paquete dplyr que facilita la realización de algunas de estas operaciones.
Partimos de una tabla de datos sintética…
df <- data.frame(1:5, letters[1:5], c(rep("a", 3), rep("b", 2)))
df## X1.5 letters.1.5. c.rep..a...3...rep..b...2..
## 1 1 a a
## 2 2 b a
## 3 3 c a
## 4 4 d b
## 5 5 e b
colnames(df) <- c("var1", "var2", "var3")
rownames(df) <- paste("subject00", 1:5, sep = "")
df## var1 var2 var3
## subject001 1 a a
## subject002 2 b a
## subject003 3 c a
## subject004 4 d b
## subject005 5 e b
df2 <- cbind(df, rnorm(5)) # añadir un vector al data frame
df2$var5 <- 5:1 # assignando valores a una nueva variable
df2## var1 var2 var3 rnorm(5) var5
## subject001 1 a a -1.304938 5
## subject002 2 b a 1.708467 4
## subject003 3 c a -1.855023 3
## subject004 4 d b 1.526245 2
## subject005 5 e b -1.995378 1
df2 <- rbind(df, list(6, "e", "b"))
df2## var1 var2 var3
## subject001 1 a a
## subject002 2 b a
## subject003 3 c a
## subject004 4 d b
## subject005 5 e b
## 6 6 e b
Existen tres formas básicas de seleccionar filas o columnas de una tablas de datos
df[1:3,]## var1 var2 var3
## subject001 1 a a
## subject002 2 b a
## subject003 3 c a
df[,c(1,3)]## var1 var3
## subject001 1 a
## subject002 2 a
## subject003 3 a
## subject004 4 b
## subject005 5 b
df[-3,-2]## var1 var3
## subject001 1 a
## subject002 2 a
## subject004 4 b
## subject005 5 b
df[,"var2"]## [1] a b c d e
## Levels: a b c d e
df$var3## [1] a a a b b
## Levels: a b
df[,c("var2","var3")]## var2 var3
## subject001 a a
## subject002 b a
## subject003 c a
## subject004 d b
## subject005 e b
df[c(T,T,F,T,F), c(T,F,T)]## var1 var3
## subject001 1 a
## subject002 2 a
## subject004 4 b
df[df[,1] == 3 | df[,3] == "b",]## var1 var2 var3
## subject003 3 c a
## subject004 4 d b
## subject005 5 e b
Un mismo conjunto de datos puede representar en una tabla de datos de acuerdo con distintas oragnizaciones o formatos.
Se denomina formato ordenado, tidy, aquella organización en la que las filas representan individuos, las columnas, variables de medidas, y las intersecciones los valores de dichas medidas.
Cuando existen más de una variable que corresponde al mismo tipo de medida, entonces los datos pueden organizarse de dos formas
El formato más adecuado dependerá de los análisis y visualizaciones que quieran hacerse.
Partiremos de un subconjunto de flights…
if(!require("nycflights13")) {
install.packages("nycflights13")
library("nycflights13")
}fl_ny2ws_W <- flights[flights$dest %in% c("IAD","BWI"),
c("origin","dest","carrier","arr_delay","dep_delay")]
fl_ny2ws_W <- cbind(key = 1:nrow(fl_ny2ws_W), fl_ny2ws_W)
head(fl_ny2ws_W)## key origin dest carrier arr_delay dep_delay
## 1 1 LGA IAD EV -14 -3
## 2 2 LGA BWI WN -19 -1
## 3 3 EWR IAD EV 12 24
## 4 4 JFK BWI MQ 851 853
## 5 5 JFK IAD B6 14 15
## 6 6 LGA IAD EV 5 10
fl_ny2ws_L <- rbind(
cbind(edge = rep("origin", nrow(fl_ny2ws_W)), fl_ny2ws_W[,c(1,4)],
airport = fl_ny2ws_W[,2], delay = fl_ny2ws_W[,6]),
cbind(edge = rep("dest", nrow(fl_ny2ws_W)), fl_ny2ws_W[,c(1,4)],
airport = fl_ny2ws_W[,3], delay = fl_ny2ws_W[,5]))
colnames(fl_ny2ws_L) <- c("edge","key","carrier","airport","delay")head(fl_ny2ws_L)## edge key carrier airport delay
## 1 origin 1 EV LGA -3
## 2 origin 2 WN LGA -1
## 3 origin 3 EV EWR 24
## 4 origin 4 MQ JFK 853
## 5 origin 5 B6 JFK 15
## 6 origin 6 EV LGA 10
tail(fl_ny2ws_L)## edge key carrier airport delay
## 14957 dest 7476 B6 IAD -4
## 14958 dest 7477 MQ BWI -9
## 14959 dest 7478 EV IAD -16
## 14960 dest 7479 EV IAD -19
## 14961 dest 7480 9E IAD -40
## 14962 dest 7481 9E BWI 2
fl_ny2ws_W2p1 <- fl_ny2ws_L[fl_ny2ws_L$edge=="origin",]
fl_ny2ws_W2p2 <- fl_ny2ws_L[fl_ny2ws_L$edge=="dest",]
fl_ny2ws_W2p1 <- fl_ny2ws_W2p1[order(fl_ny2ws_W2p1$key),-1]
fl_ny2ws_W2p2 <- fl_ny2ws_W2p2[order(fl_ny2ws_W2p2$key),-c(1,2)]
fl_ny2ws_W2 <- cbind(fl_ny2ws_W2p1,fl_ny2ws_W2p2[,-1])
colnames(fl_ny2ws_W2) <- c("key", "carrier", "origin", "dep_delay",
"dest","arr_delay")head(fl_ny2ws_W2)## key carrier origin dep_delay dest arr_delay
## 1 1 EV LGA -3 IAD -14
## 2 2 WN LGA -1 BWI -19
## 3 3 EV EWR 24 IAD 12
## 4 4 MQ JFK 853 BWI 851
## 5 5 B6 JFK 15 IAD 14
## 6 6 EV LGA 10 IAD 5
tail(fl_ny2ws_W2)## key carrier origin dep_delay dest arr_delay
## 7476 7476 B6 JFK 12 IAD -4
## 7477 7477 MQ JFK -5 BWI -9
## 7478 7478 EV LGA -5 IAD -16
## 7479 7479 EV JFK -1 IAD -19
## 7480 7480 9E JFK -7 IAD -40
## 7481 7481 9E JFK 38 BWI 2
La forma más habitual de eliminar filas o columnas es segmentando la tabla de datos. Sin embargo, una columna también puede eliminarse asignando la misma a NULL.
df <- data.frame(1:5, letters[1:5], c(rep("a", 3), rep("b", 2)))
colnames(df) <- c("var1", "var2", "var3")
rownames(df) <- paste("subject00", 1:5, sep = "")
df <- df[-2,]
df$var2 <- NULLdf## var1 var3
## subject001 1 a
## subject003 3 a
## subject004 4 b
## subject005 5 b
Si lo que se desea es eliminar una variable del entorno de trabajo entonces se usa la función rm.
rm(df)sum_fl_ny2ws <- data.frame(edge=c("origin","dest"))
sum_fl_ny2ws$mean_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
mean,na.rm=TRUE)
sum_fl_ny2ws$r_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
function(x) diff(range(x,na.rm=TRUE)))
sum_fl_ny2ws$s_delay <- apply(fl_ny2ws_W[,c("dep_delay","arr_delay")],2,
sd,na.rm=TRUE)
sum_fl_ny2ws## edge mean_delay r_delay s_delay
## 1 origin 16.84267 885 47.85548
## 2 dest 13.11556 904 51.52259
sum_fl_ny2ws <- data.frame(edge=c("origin","dest"))
sum_fl_ny2ws$mean_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,
mean,na.rm=TRUE)
sum_fl_ny2ws$r_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,
function(x) diff(range(x,na.rm=TRUE)))
sum_fl_ny2ws$s_delay <- by(fl_ny2ws_L$delay,fl_ny2ws_L$edge,sd,na.rm=TRUE)
sum_fl_ny2ws## edge mean_delay r_delay s_delay
## 1 origin 16.84267 885 47.85548
## 2 dest 13.11556 904 51.52259
dplyrEstas mismas operaciones que se han comentado más arriba, se pueden llevar a cabo también a partir de un sistema coherente de métodos creado como una gramática para la manipulación de datos, el paquete dplyr –que forma parte de tidyverse–.
En dplyr las operaciones básicas, se realizan mediante cuatro métodos:
select: segmentar columnas,filter: segmentar filas,mutate: crear nuevas columnas, yarrange: ordenar.Estas se encadenan usando el operador pipe, %>%.
df <- flights %>% select(origin, dest, arr_delay) %>%
filter(origin == "LGA" & (dest == "IAD" | dest == "BWI")) %>%
mutate(arr_delay_h=arr_delay/60) %>%
arrange(-arr_delay_h)
df## # A tibble: 1,818 x 4
## origin dest arr_delay arr_delay_h
## <chr> <chr> <dbl> <dbl>
## 1 LGA IAD 398 6.633333
## 2 LGA IAD 326 5.433333
## 3 LGA IAD 312 5.200000
## 4 LGA IAD 306 5.100000
## 5 LGA IAD 292 4.866667
## 6 LGA IAD 281 4.683333
## 7 LGA IAD 280 4.666667
## 8 LGA IAD 265 4.416667
## 9 LGA IAD 262 4.366667
## 10 LGA IAD 253 4.216667
## # ... with 1,808 more rows
Para la creación de resúmenes a partir de tablas en formato largo, es muy útil y cómoda la combinación group_by y summarize.
df %>% group_by(dest) %>% summarise(mean_delay = mean(arr_delay, na.rm=TRUE))## # A tibble: 2 x 2
## dest mean_delay
## <chr> <dbl>
## 1 BWI -7.866667
## 2 IAD 13.708258
Por último, las conversiones entre formatos también pueden hacerse a partir de las funciones spread y gather del paquete tidyr –incluido en tidyverse y que se integra de forma natural en la gramática propuestas por dplyr–.